Mathematical Biosciences

Cancer-associated fibroblasts (CAFs) are a central component of the tumor microenvironment that facilitate a supportive niche for cancer progression and metastasis. Experimental evidence suggests that CAFs can facilitate estrogen-independent tumor growth, thereby reducing the efficacy of anti-hormonal therapies. Understanding and quantifying the complex interactions between tumor cells, hormonal signalling, and the microenvironment are crucial for designing more effective and individualized treatment strategies. We propose a mathematical framework to explore the influence of CAFs on ER+ breast cancer progression and to evaluate strategies to mitigate their impact. We develop a system of nonlinear ordinary differential equations that substantiates the experimental observations by providing a mechanistic basis for the role of CAFs in regulating estrogen-independent growth dynamics. We then employ optimal control theory to evaluate distinct therapeutic approaches involving monotherapy or combinations of: (i) inhibition of tumor-to-CAF signaling, (ii) inhibition of CAF-to-tumor proliferative signaling, and (iii) endocrine therapy. Taken together, our results demonstrate that CAF-targeted strategies can enhance treatment efficacy across various estrogen dosing regimens. Our study provides new insights into the potential of CAF as a therapeutic target that could help to improve existing approaches for endocrine therapies.

Vaccine-acquired immunity plays an important role in controlling the spread of many infectious diseases; however, vaccine efficacy can diminish over time. This work uses a mathematical model to study the effects of waning vaccination-acquired immunity on infection incidence. With an SEIR-type compartmental model that considers both vaccinated and unvaccinated populations (and their mixing), we present mathematical conditions under which vaccinated individuals drive ongoing growth in infections, i.e., over half of the daily incidence arises from vaccinated individuals. Analysis of a mathematical model of COVID-19 spread in the state of Colorado suggests how and for what duration vaccinated individuals could have sustained such growth. Importantly, our model demonstrates that, despite potential for brief vaccinated-driven periods of growth in infections, which occur among unvaccinated-driven periods of growth in infections, increased vaccination coverage always reduces total cases and total hospitalizations. This work provides insight into how waning immunity in vaccinated populations can contribute to ongoing infection incidence and demonstrates the value of complementary interventions to prevent disease spread in vaccinated populations.

Evolutionary therapies regulate heterogeneous populations by altering selective pressures through treatment sequences in cancer and infections. This letter develops an invariant-set framework for treatment-induced containment based on positive triangular invariant sets. For periodically switched systems, sufficient conditions are derived for the existence of such invariant regions. Robustness with respect to mutation is established by showing that the invariant simplex persists under small perturbations of the subsystem matrices. In the two-phenotype case, the analysis yields an explicit mutation threshold that separates regimes in which therapy cycling maintains containment from regimes in which mutation can enable evolutionary escape. Simulations illustrate the geometry of the invariant sets and the role of mutation and dwell time in containment robustness.

Canine rabies persists in NDjamena (Chad) despite vaccination campaigns exceeding 70% coverage, suggesting a role for dog mobility and spatial heterogeneity. We propose a metapopulation SEIR model incorporating distance-modulated dog movements and an explicit vaccinated class. Analysis of the isolated patch establishes global stability of the disease-free equilibrium via a Lyapunov function. For the metapopulation, a composite Lyapunov function shows that elimination is governed by a reproduction number [R]v. Calibrated with field data (2012-2022), simulations reveal that uniform vaccination of both patches reduces [R]v by 46% (from 2.84 to 1.52) but does not achieve elimination, while targeted strategies are less effective. These results demonstrate that exhaustive vaccination coverage across the entire urban network and increased vaccination intensity are necessary to eliminate canine rabies in NDjamena. Our model provides a quantitative framework for planning effective control strategies.

Coinfection of COVID-19 and malaria in endemic regions may generate complex epidemiological interactions that influence susceptibility patterns, disease burden, and outbreak risk. Although malaria-acquired immunity has been hypothesized to modulate host responses to other infections, its population-level implications for COVID-19 transmission under uncertainty remain insufficiently understood. In this study, we develop a deterministic-stochastic compartmental model for the coupled dynamics of COVID-19, malaria, and their co-infection. Malaria-acquired partial immunity is incorporated through a relative susceptibility parameter that reduces the risk of COVID-19 infection among malaria-recovered individuals. For the deterministic system, we establish positivity, boundedness, an invariant feasible region, and basic reproduction numbers for the COVID-19-only and malaria-only subsystems. We then use numerical simulations to examine how immunity-mediated reductions in susceptibility may influence COVID-19 incidence, peak burden, hospitalization, and cumulative mortality. To account for environmental and transmission variability, we extend the deterministic model to an Ito stochastic differential equation framework and use repeated realizations to characterize uncertainty in epidemic trajectories, peak distributions, and outbreak risk. In addition, global sensitivity analysis based on partial rank correlation coefficients (PRCCs) is performed to identify the parameters with the greatest influence on COVID-19 outcomes. Our results suggest that, under the assumed modeling framework, malaria-acquired partial immunity may reduce the peak infectious burden and cumulative mortality associated with COVID-19. The stochastic simulations further show substantial variability around deterministic trajectories and indicate a non-negligible probability of large outbreak events that are not fully captured by mean-field predictions alone. Overall, the proposed framework provides an uncertainty-aware, mechanistic basis for studying COVID-19-malaria co-dynamics and for assessing how interacting disease processes may shape epidemic outcomes in endemic settings.

Human immunodeficiency virus (HIV) persistence remains a major barrier to cure due to the existence of long-lived latent reservoirs that evade immune clearance and persist despite combination antiretroviral therapy (ART). Although ART effectively suppresses viral replication, treatment interruption often leads to rapid viral rebound originating from these latent reservoirs. In this study, we develop a deterministic mathematical model describing the in vivo dynamics of HIV infection incorporating uninfected CD4+ T cells, infected cells, latent reservoirs, deep latent reservoirs, and infectious and non-infectious virions, while explicitly accounting for the therapeutic effects of reverse transcriptase inhibitors (RTIs), protease inhibitors (PIs), and Tat transcription inhibitors. Analytical results establish positivity and boundedness of solutions and derive the effective reproduction number Re using the next-generation matrix approach. Stability analysis shows that the virus-free equilibrium is locally asymptotically stable when Re < 1, while viral persistence occurs when Re > 1. Numerical simulations were performed to investigate therapy interactions, viral rebound following treatment interruption, and the impact of drug efficacy on viral set-points and latent reservoir dynamics. To further explore therapy interactions, three-dimensional viral set-point surfaces and heat maps were generated to examine how combinations of infection inhibition, viral production inhibition, and transcriptional inhibition influence viral dynamics. The simulations reveal that Tat inhibition suppresses viral transcription, thereby reducing the transition of infected cells into productive infection and limiting viral propagation when combined with conventional ART mechanisms. The therapy parameter planes further demonstrate that strong transcriptional inhibition promotes the transition of infected cells into deep latency, supporting the emerging block-and-lock strategy for functional HIV cure. In addition, a three-dimensional eradication boundary surface and therapy cube were constructed to identify regions of parameter space where Re < 1, corresponding to successful viral control. These visualizations show that viral eradication is unlikely when therapies act independently but becomes achievable when multiple therapeutic mechanisms act simultaneously. Overall, the results highlight the critical role of transcriptional inhibition through Tat-targeting therapies in complementing existing ART regimens. By simultaneously suppressing viral replication and promoting deep latency, Tat-based combination strategies may significantly reduce viral rebound and contribute to long-term functional control of HIV infection.

AO_SCPLOWBSTRACTC_SCPLOWOncolytic viruses provide cancer therapy using replication-competent viruses that selectively infect and lyse tumour cells. Their tumour-specific replication also enables the delivery of targeted, virus-encoded gene products, such as enzymes that activate prodrugs. This dual functionality offers the potential for synergistic effects by combining direct oncolysis with localised drug activation. The interplay between infection, replication, lysis, and gene product delivery remains poorly understood. Here, we introduce a spatially structured, multi-patch model of cancer cells infected by an oncolytic virus engineered to deliver a prodrug-activating enzyme. The spatial system is first represented as a microscopic model and subsequently reduced via spectral dimension reduction techniques. This reduction yields an ordinary differential equation model for a set of coarse-grained variables, which we analyze both without the transgene (OV model) and with the transgene (TOV model). For each case, we compute the basic reproduction number, R0, which determines the conditions for viral spread. Both models exhibit three regimes via transcritical bifurcations: (i) R0 < 0, extinction of both cancer and infected cells; (ii) 0 < R0 [&le;] 1, persistence of cancer cells only; and (iii) R0 > 1, coexistence as a stable node or as a focus. The TOV model, as a difference form the OV model, can undergo periodic oscillations arising from a Hopf-Andronov bifurcation. Notably, oscillation amplitudes can be controlled such that cancer cells largely decrease when drug-induced death is stronger in non-infected cells than in infected ones, enabling effective cancer cells killing while maintaining viral replication and prodrug activation. The qualitative behaviour of the coarse-grained model is shown to be preserved in both the microscopic and spatially explicit models.

Cells communicate via extracellular ligands, such as hormones, which bind to plasma membrane receptors and trigger intracellular signaling cascades. G Protein-Coupled Receptors (GPCRs) exemplify this mechanism by initiating signaling both at the cell surface and, from intracellular compartments such as endosomes. The kinetics and spatial localization of these signals are critical determinants of cellular responses, yet receptor trafficking-including internalization, endosomal sorting, and recycling-remains a pivotal but often overlooked component of theoretical GPCR models. In this study, we present a mathematical framework that integrates receptor trafficking and signaling compartmentalization into generic GPCR dynamic models. Using a compartmentalized approach based on systems of ordinary differential equations (Chemical Reaction Networks), we analyze how receptor internalization and recycling modulate ligand-induced responses. Our results show that the balance between plasma membrane and endosomal signaling can significantly enhance or diminish ligand efficacy. Calibrated with high-throughput kinetic data, our model offers a refined tool for ligand pharmacological characterization and advances the understanding of GPCR signaling spatial organization.

Successful development of multicellular organisms requires cells to transition between states with precise timing. Distinct cell states are often understood as being maintained by stabilizing regulatory networks, such that a complete cell-state transition requires network rewiring through partial dismantling of the current state and concurrent reconfiguration into a new one. Empirically, these transitions are often investigated by quantifying the gain or loss of expression of a small number of state-specific markers, frequently a single proxy. A general quantitative framework for comparing the kinetics of such transitions across experimental conditions is lacking. Here, we show that the delayed Weibull distribution provides a natural description of cell-state transition kinetics when transition is viewed as the cumulative consequence of many molecular events, whose timing may vary between cells and conditions, analogous to system failure in reliability theory. This formulation yields a compact model with interpretable parameters describing the delay before transition onset, the characteristic timescale of transition, the temporal form of the transition hazard, and the fraction of cells competent to respond. Together, this framework provides a practical and interpretable approach for quantifying the kinetics of cell-state transitions and how they are altered by perturbation.

Variability in gene expression among single cells and growing cell populations can arise from the stochastic nature of protein synthesis, which often occurs in random bursts. This study investigates the variability in the expression of a growth-sustaining protein, whose concentration is regulated by a negative feedback loop due to cell growth-induced dilution. We model the distribution of protein concentration using a Chapman-Kolmogorov equation for single cells and a population balance equation for growing cell populations. For single cells, we derive an explicit solution for the protein concentration distribution in state space and represent it as a Bessel function in Laplace space. For growing populations, we find that the distribution satisfies a Heun differential equation with singular boundary conditions. By addressing the central connection problem for the Heun equation, we quantify the population-level protein distribution and determine the Mathusian parameter, which characterizes population growth. This work provides a comprehensive analytical framework for understanding how stochastic protein synthesis impacts gene expression variability and population dynamics.

Metabolic networks are typically viewed as homeostatic systems that stabilize flux, energy charge, redox balance, and metabolite availability under perturbation. However, it remains unclear whether the same feedback architectures that support metabolic robustness can also generate learning-like, experience-dependent adaptation. Here, we develop a coarse-grained dynamical model of mammalian energy metabolism to test whether prior perturbation can improve future metabolic responses. The model represents core glucose, glutamine, fatty acid, and oxidative phosphorylation pathways as coupled ordinary differential equations with Michaelis-Menten-type fluxes, product-inhibition feedback, adaptive enzyme-capacity regulation, and explicit ATP costs for enzyme adjustment. Rather than aiming to reproduce quantitative fluxes for a specific cell type, the framework is designed to expose how metabolic feedback, regulatory cost, repeated perturbation, and environmental variability interact. We use this model to ask whether adaptive enzyme regulation enables improved recovery after repeated challenges, whether such effects depend on energetic control costs, and whether environmental variability broadens or constrains the set of reachable adaptive states. This approach provides a tractable way to investigate how homeostatic metabolic regulation may give rise to experience-dependent metabolic plasticity.

Ordinary differential equation models of biochemical reactions are often formulated as stoichiometric systems in which the dynamics arise from a collection of interacting processes. A central challenge is that the functional form of each process is rarely known a priori and may be difficult to infer from data. We propose biochemically informed neural ordinary differential equations (BINODEs), a neural-ODE framework that retains the stoichiometric structure of mechanistic models while representing individual processes by neural networks. In BINODEs, the outputs of neural network processes (NNPs) are mapped to state derivatives through a linear layer analogous to a stoichiometric matrix. This architecture allows biological side information, such as process-specific inputs, sign constraints, and monotonicity assumptions, to be built directly into the model. We characterize the approximation properties of NNPs for several standard biochemical rate laws and show that the proposed framework recovers both trajectories and process-level structure in Monod, Lotka-Volterra, pharmacokinetic, and ultradian endocrine models. These results suggest that BINODEs offer a useful compromise between mechanistic interpretability and data-driven flexibility for modeling partially known biochemical or biological dynamical systems.

A novel mathematical framework to define the threshold of action potentials in excitable cells is presented. Unlike previously applied methods that rely on approximations or specific fixed-point bifurcations, the approach focuses on the geometry of membrane potential trajectories. Specifically, the focus is on the concavity changes during the upstroke of an electrical pulse. These changes in concavity form a curve of inflection points that defines a region in phase space crossed by all the action potentials in the system, and containing no non-action potential trajectories. Such region is called the excitability region and its size can be measured, thus providing a measure for the excitability of a dynamical system, and a way to compare the excitability between systems representing different biological phenotypes and stimulus conditions. The work transforms the traditionally vague physiological concept of excitability into a rigorous analytical description applicable across continuous, single compartment models of electrical excitability.

This work focuses on the global and partial identification problem for fractional differential equations. We provide a general numerical procedure based on global and local optimization algorithms with two refinements for biological systems that ensure solution positivity and homogeneous parameter units. The method is applied to a new fractional model of Dengue outbreak called the Fractional Homogeneous Nishiura (FHN) model, calibrated using data of newly infected people in Cape Verde. We show that our identification method yields a better fit between data and model solutions than previous approaches and that our FHN model captures the dynamics of Dengue more closely than existing systems.

Intra-tumor heterogeneity in proliferation rates fundamentally influences cancer progression and treatment resistance. To investigate how continuous phenotypic variation shapes eco-evolutionary dynamics, we develop a phenotype-structured partial differential equation framework that explicitly models proliferation het-erogeneity as a dynamic trait distribution. Our model integrates three key biological principles: (1) phenotypic diffusion capturing heritable variation in proliferation rates, (2) global resource competition enforcing density-dependent growth constraints, and (3) an experimentally grounded life-history trade-off linking elevated proliferation to increased mortality. Using adaptive dynamics, we derive the optimum proliferation rate in a growing tumor, showing that the optimal phenotype dynamically shifts toward slower proliferation as tumors approach carrying capacity. We perform in silico treatment simulations for four different treatment regimes (pan-proliferation, low-, mid-, and high-proliferation targeting) to show how therapeutic selective pressures reshape fitness landscapes. While all treatments slow down tumor growth, they induce divergent evolutionary trajectories: low- and mid-proliferation targeting enrich fast-proliferating clones, whereas high-proliferation targeting selects for slower phenotypes. We connect these dynamics with changes in mean proliferation rates during and after treatment. We use adaptive dynamics to explain the shifts in mean proliferation rate during treatment, showing how each regimen alters the maximum fitness proliferation rate. Our work establishes a predictive, evolutionarily grounded framework for understanding how therapy reshapes tumor proliferation landscapes, offering a mechanistic basis for designing strategies that anticipate and counteract adaptive resistance.

Tetanus is a severe disease of the nervous system, transmitted through bacteria in the environment. In the absence of medical attention, case fatality rates are extremely high. Despite progress towards maternal and neonatal tetanus elimination targets, tetanus remains a serious public health problem. Routine infant and maternal vaccination have contributed to considerable reduction in cases and deaths from tetanus globally. However, protective immunity wanes over time. To increase duration of protection, the World Health Organization recommends three diphtheria-tetanus-pertussis-containing vaccine booster doses be given in early childhood, childhood, and adolescence. Evidence to support country-level decision-making about the introduction of these booster doses is critical. We have developed a novel age-structured, deterministic compartmental model of tetanus transmission and vaccination. The model is driven by environmental transmission and incorporates interventions like hygiene and safe birth practices to reduce the magnitude of environmental transmission. It explicitly models vaccination, separating each dose of the primary series, booster doses, and maternal vaccination to capture dose-specific effectiveness and duration of protection. The model captures heterogeneous immunity profiles by dose and age, and the cumulative nature of vaccine-derived protection. The immune dynamics follow the patterns described in literature and can replicate seroprevalence studies, although the exact characterisation of immunity in the literature still has gaps. This model presents a substantial advancement on previously published models and is well positioned to inform tailored vaccination strategies to reduce neonatal and non-neonatal tetanus.

Parameter estimation in nonlinear biological dynamical systems is a difficult inverse problem because the governing equations are often stiff or oscillatory, the data are sparse and noisy, and the objective landscape is non-convex. Physics-informed neural networks (PINNs) offer an alternative to purely simulation-based calibration by representing state trajectories with neural networks while penalizing violations of the governing equations. This paper studies the empirical reliability of PINNs for recovering the parameters of the repressilator, a synthetic genetic oscillator formed by three cyclically repressive genes. We use synthetic time-series generated from the standard ordinary differential equation model and train inverse PINNs to estimate the production parameter {beta} and the Hill coefficient n. The study varies observation noise, partial observation of repressors, sampling density, sensitivity to initial parameter guesses, and the difference between stable and oscillatory regimes. The results show that PINNs can reconstruct trajectories accurately when the model structure is correct and the three repressors are observed, but parameter recovery is more fragile than trajectory fitting. Noise, sparse sampling, unobserved variables, and unfavorable initial guesses increase the risk of biased estimates. The stable regime is easier to reconstruct, whereas the oscillatory regime provides richer information but also exposes optimization sensitivity. These findings support PINNs as a useful reverse-engineering tool for small gene-regulatory ODE models, while highlighting the need for repeated runs, uncertainty reporting, and experimental designs that improve identifiability.

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

Mass-conserved reaction-diffusion systems are used as mathematical models for various phenomena such as cell polarity. Numerical simulations of this system present transient dynamics in which multiple stripe patterns converge to spatially monotonic patterns. Previous studies indicated that the transient dynamics are driven by a mass conservation law and by variations in the amount of substance contained in each pattern, which we refer to as "pattern flux". However, it is challenging to mathematically investigate these pattern dynamics. In this study, we introduce a reaction-diffusion compartment model to investigate the pattern dynamics in view of the conservation law and the pattern flux. This model is defined on multiple intervals (compartments), and diffusive couplings are imposed on each boundary of the compartments. Corresponding to the transient dynamics in the original system, we consider the dynamics around stripe patterns in the compartment model. We derive ordinary differential equations describing the pattern dynamics of the compartment model and analyze the existence and stability of equilibria for the reduced ODE with respect to the boundary parameters. For a specific parameter setting, we obtained results consistent with previous studies. Moreover, we present that the stripe patterns in the compartment model are potentially stabilized by changing the parameter, which is not observed in the original system. We expect that the methodology developed in this paper is extendable to various directions, such as membrane-induced pattern control.

Malaria persists worldwide, exerting its greatest impact in sub-Saharan Africa. This study develops and uses a mathematical model to assess how sub-optimum versus optimum treatment of malaria drives asymptomatic infections, immunity build-up, and sustained transmission, providing insights for effective control. Fitted to case data from Kenya and Nigeria, the framework is used to quantify the burden of malaria and the additional cost associated with sub-optimum treatment. Global sensitivity analysis identifies mosquito demographic parameters, biting rates, and malaria treatment rate among major disease drivers under sub-optimum treatment, emphasizing the need for integrated strategies that improve access to optimum treatment and reduce mosquito-human contact. Model simulations show that sub-optimum treatment amplifies asymptomatic prevalence, sustaining/increasing malaria transmission and burden. Further simulations reveal that optimum treatment could avert more than one-third of infections and deaths, while asymptomatic infections contribute up to $96%$ ($75%$) of malaria-related Years Lived with Disability in Kenya (Nigeria). Cost analysis shows that optimum treatment lowers malaria burden significantly and can reduce annual total treatment costs by $\approx $12$ million, underscoring the substantial economic and public health gains of limiting sub-optimum care. This study demonstrates that effective and sustained malaria control requires strengthening adherence to treatment, minimizing sub-optimum treatment, reducing mosquito-human contact, and targeting asymptomatic carriers to curb hidden transmission and reduce long-term health and economic losses.